#include "luminance.wgsl"
#include "space/srgb2rgb.wgsl"
#include "space/rgb2srgb.wgsl"
#include "space/xyz2srgb.wgsl"

fn mixSpectral(srgb1: vec3f, srgb2: vec3f, t: f32) -> vec3f {
    let lrgb1 = srgb1;
    let lrgb2 = srgb2;

    // Linear luminance to concentration
    let t1 = luminanceLinear(lrgb1) * pow(1.0 - t, 2.0);
    let t2 = luminanceLinear(lrgb2) * pow(t, 2.0);
    let pct = t2 / (t1 + t2);

    var R1: array<f32,37>;
    R1[0] = dot(vec3f(0.03065266, 0.00488428, 0.96446343), lrgb1);
    R1[1] = dot(vec3f(0.03065266, 0.00488428, 0.96446661), lrgb1);
    R1[2] = dot(vec3f(0.03012503, 0.00489302, 0.96499804), lrgb1);
    R1[3] = dot(vec3f(0.02837440, 0.00505932, 0.96660443), lrgb1);
    R1[4] = dot(vec3f(0.02443079, 0.00552416, 0.97009698), lrgb1);
    R1[5] = dot(vec3f(0.01900359, 0.00668451, 0.97438902), lrgb1);
    R1[6] = dot(vec3f(0.01345743, 0.00966823, 0.97695626), lrgb1);
    R1[7] = dot(vec3f(0.00905147, 0.01843871, 0.97257598), lrgb1);
    R1[8] = dot(vec3f(0.00606943, 0.05369084, 0.94029002), lrgb1);
    R1[9] = dot(vec3f(0.00419240, 0.30997719, 0.68585683), lrgb1);
    R1[10] = dot(vec3f(0.00300621, 0.84166297, 0.15533920), lrgb1);
    R1[11] = dot(vec3f(0.00229452, 0.95140393, 0.04629964), lrgb1);
    R1[12] = dot(vec3f(0.00190474, 0.97711658, 0.02096763), lrgb1);
    R1[13] = dot(vec3f(0.00175435, 0.98538119, 0.01284767), lrgb1);
    R1[14] = dot(vec3f(0.00182349, 0.98819579, 0.00996131), lrgb1);
    R1[15] = dot(vec3f(0.00218287, 0.98842729, 0.00937163), lrgb1);
    R1[16] = dot(vec3f(0.00308472, 0.98651266, 0.01038752), lrgb1);
    R1[17] = dot(vec3f(0.00539517, 0.98125477, 0.01334010), lrgb1);
    R1[18] = dot(vec3f(0.01275154, 0.96796653, 0.01927821), lrgb1);
    R1[19] = dot(vec3f(0.04939664, 0.92126320, 0.02934328), lrgb1);
    R1[20] = dot(vec3f(0.41424516, 0.54678757, 0.03897609), lrgb1);
    R1[21] = dot(vec3f(0.89425217, 0.07097922, 0.03478168), lrgb1);
    R1[22] = dot(vec3f(0.95202201, 0.02151275, 0.02647979), lrgb1);
    R1[23] = dot(vec3f(0.96833286, 0.01120932, 0.02047109), lrgb1);
    R1[24] = dot(vec3f(0.97175685, 0.00778212, 0.02047109), lrgb1);
    R1[25] = dot(vec3f(0.97320302, 0.00633303, 0.02047109), lrgb1);
    R1[26] = dot(vec3f(0.97387285, 0.00566048, 0.02047109), lrgb1);
    R1[27] = dot(vec3f(0.97418395, 0.00534751, 0.02047109), lrgb1);
    R1[28] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[29] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[30] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[31] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[32] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[33] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[34] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[35] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);
    R1[36] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb1);

    var R2: array<f32,37>;
    R2[0] = dot(vec3f(0.03065266, 0.00488428, 0.96446343), lrgb2);
    R2[1] = dot(vec3f(0.03065266, 0.00488428, 0.96446661), lrgb2);
    R2[2] = dot(vec3f(0.03012503, 0.00489302, 0.96499804), lrgb2);
    R2[3] = dot(vec3f(0.02837440, 0.00505932, 0.96660443), lrgb2);
    R2[4] = dot(vec3f(0.02443079, 0.00552416, 0.97009698), lrgb2);
    R2[5] = dot(vec3f(0.01900359, 0.00668451, 0.97438902), lrgb2);
    R2[6] = dot(vec3f(0.01345743, 0.00966823, 0.97695626), lrgb2);
    R2[7] = dot(vec3f(0.00905147, 0.01843871, 0.97257598), lrgb2);
    R2[8] = dot(vec3f(0.00606943, 0.05369084, 0.94029002), lrgb2);
    R2[9] = dot(vec3f(0.00419240, 0.30997719, 0.68585683), lrgb2);
    R2[10] = dot(vec3f(0.00300621, 0.84166297, 0.15533920), lrgb2);
    R2[11] = dot(vec3f(0.00229452, 0.95140393, 0.04629964), lrgb2);
    R2[12] = dot(vec3f(0.00190474, 0.97711658, 0.02096763), lrgb2);
    R2[13] = dot(vec3f(0.00175435, 0.98538119, 0.01284767), lrgb2);
    R2[14] = dot(vec3f(0.00182349, 0.98819579, 0.00996131), lrgb2);
    R2[15] = dot(vec3f(0.00218287, 0.98842729, 0.00937163), lrgb2);
    R2[16] = dot(vec3f(0.00308472, 0.98651266, 0.01038752), lrgb2);
    R2[17] = dot(vec3f(0.00539517, 0.98125477, 0.01334010), lrgb2);
    R2[18] = dot(vec3f(0.01275154, 0.96796653, 0.01927821), lrgb2);
    R2[19] = dot(vec3f(0.04939664, 0.92126320, 0.02934328), lrgb2);
    R2[20] = dot(vec3f(0.41424516, 0.54678757, 0.03897609), lrgb2);
    R2[21] = dot(vec3f(0.89425217, 0.07097922, 0.03478168), lrgb2);
    R2[22] = dot(vec3f(0.95202201, 0.02151275, 0.02647979), lrgb2);
    R2[23] = dot(vec3f(0.96833286, 0.01120932, 0.02047109), lrgb2);
    R2[24] = dot(vec3f(0.97175685, 0.00778212, 0.02047109), lrgb2);
    R2[25] = dot(vec3f(0.97320302, 0.00633303, 0.02047109), lrgb2);
    R2[26] = dot(vec3f(0.97387285, 0.00566048, 0.02047109), lrgb2);
    R2[27] = dot(vec3f(0.97418395, 0.00534751, 0.02047109), lrgb2);
    R2[28] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[29] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[30] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[31] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[32] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[33] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[34] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[35] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);
    R2[36] = dot(vec3f(0.97432335, 0.00520568, 0.02047109), lrgb2);

    var R: array<f32,37>;
    for (var i = 0; i < 37; i++) {
      var KS = 0.0;
      KS += (1.0 - pct) * (pow(1.0 - R1[i], 2.0) / (2.0 * R1[i]));
      KS += pct * (pow(1.0 - R2[i], 2.0) / (2.0 * R2[i]));
      let KM = 1.0 + KS - sqrt(pow(KS, 2.0) + 2.0 * KS);
      //Saunderson correction
      R[i] = KM;
    }

    var xyz = vec3f(0.0);
    xyz +=  R[0] * vec3f(0.00013656, 0.00001886, 0.00060998);
    xyz +=  R[1] * vec3f(0.00131637, 0.00018140, 0.00595953);
    xyz +=  R[2] * vec3f(0.00640948, 0.00080632, 0.02967913);
    xyz +=  R[3] * vec3f(0.01643026, 0.00203723, 0.07855475);
    xyz +=  R[4] * vec3f(0.02407799, 0.00344701, 0.11921905);
    xyz +=  R[5] * vec3f(0.03573573, 0.00662872, 0.18425721);
    xyz +=  R[6] * vec3f(0.03894236, 0.01029015, 0.21077492);
    xyz +=  R[7] * vec3f(0.03004572, 0.01410577, 0.17634280);
    xyz +=  R[8] * vec3f(0.01860940, 0.01944240, 0.12741269);
    xyz +=  R[9] * vec3f(0.00745020, 0.02631783, 0.07374910);
    xyz += R[10] * vec3f(0.00129006, 0.03273183, 0.03663768);
    xyz += R[11] * vec3f(0.00052314, 0.04424704, 0.01923495);
    xyz += R[12] * vec3f(0.00344737, 0.05701200, 0.00807728);
    xyz += R[13] * vec3f(0.01065677, 0.06907721, 0.00335702);
    xyz += R[14] * vec3f(0.02169564, 0.08047999, 0.00140323);
    xyz += R[15] * vec3f(0.03395004, 0.08541136, 0.00053757);
    xyz += R[16] * vec3f(0.04732762, 0.08725039, 0.00020463);
    xyz += R[17] * vec3f(0.06029657, 0.08416902, 0.00007431);
    xyz += R[18] * vec3f(0.07284094, 0.07860677, 0.00002725);
    xyz += R[19] * vec3f(0.08385845, 0.07114656, 0.00001053);
    xyz += R[20] * vec3f(0.08612109, 0.05907490, 0.00000391);
    xyz += R[21] * vec3f(0.08746894, 0.05050107, 0.00000166);
    xyz += R[22] * vec3f(0.07951403, 0.04005938, 0.00000072);
    xyz += R[23] * vec3f(0.06405614, 0.02932589, 0.00000000);
    xyz += R[24] * vec3f(0.04521591, 0.01939909, 0.00000000);
    xyz += R[25] * vec3f(0.03062648, 0.01258803, 0.00000000);
    xyz += R[26] * vec3f(0.01838938, 0.00734245, 0.00000000);
    xyz += R[27] * vec3f(0.01044320, 0.00409671, 0.00000000);
    xyz += R[28] * vec3f(0.00576692, 0.00223674, 0.00000000);
    xyz += R[29] * vec3f(0.00279715, 0.00107927, 0.00000000);
    xyz += R[30] * vec3f(0.00119535, 0.00046015, 0.00000000);
    xyz += R[31] * vec3f(0.00059496, 0.00022887, 0.00000000);
    xyz += R[32] * vec3f(0.00029365, 0.00011300, 0.00000000);
    xyz += R[33] * vec3f(0.00011500, 0.00004432, 0.00000000);
    xyz += R[34] * vec3f(0.00006279, 0.00002425, 0.00000000);
    xyz += R[35] * vec3f(0.00003275, 0.00001268, 0.00000000);
    xyz += R[36] * vec3f(0.00001376, 0.00000535, 0.00000000);

    return srgb2rgb(xyz2srgb(xyz));
}